rm(list = ls())
require(foreign)
require(stats4)
require(extRemes)
require(stargazer)

data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)

dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)


#some descriptive stats and graphs
exp2frame=rbind(dataframe21,dataframe22,dataframe23)
exp2frame=subset(exp2frame,uid>800)
exp2frame$correct <-(exp2frame$state==4&(exp2frame$response==2|exp2frame$response==7|exp2frame$response==8|exp2frame$response==9))|(exp2frame$state==3&(exp2frame$response==1|exp2frame$response==6|exp2frame$response==3|exp2frame$response==5))
exp2frame$incentive <- (exp2frame$qid==5)*5+(exp2frame$qid==1)*40+(exp2frame$qid==6)*70+(exp2frame$qid==7)*95
accuracyexp2=c(sum((exp2frame$incentive==5)*exp2frame$correct)/sum(exp2frame$incentive==5),sum((exp2frame$incentive==40)*exp2frame$correct)/sum(exp2frame$incentive==40),sum((exp2frame$incentive==70)*exp2frame$correct)/sum(exp2frame$incentive==70),sum((exp2frame$incentive==95)*exp2frame$correct)/sum(exp2frame$incentive==95))
incentivec=c(5,40,70,95)
uidlistexp2=unique(exp2frame$uid)

#check N
#**
length(uidlistexp2)
#**

exp2totcor=c(sum(exp2frame$correct*(exp2frame$incentive==5)),sum(exp2frame$correct*(exp2frame$incentive==40)),sum(exp2frame$correct*(exp2frame$incentive==70)),sum(exp2frame$correct*(exp2frame$incentive==95)))
exp2totinc=c(sum((1-exp2frame$correct)*(exp2frame$incentive==5)),sum((1-exp2frame$correct)*(exp2frame$incentive==40)),sum((1-exp2frame$correct)*(exp2frame$incentive==70)),sum((1-exp2frame$correct)*(exp2frame$incentive==95)))


#find the multinomial constant
multinom=function(vec){
if(length(vec)==2){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
}
if(length(vec)==3){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2]),1:max(1,vec[3])))))
}
}



exp2multinomvec=vector("numeric")
for(i in 1:4){
exp2multinomvec[i]=multinom(c(exp2totcor[i],exp2totinc[i]))
}
exp2multinomconst=sum((exp2multinomvec))

######################
#linear mutual information

exp2acculin=function(kappaall, kappaloc){
accuvec=vector('numeric')
for(i in 1:4){
cost=function(accu){
return((kappaall+kappaloc)*(accu*log(accu)+(1-accu)*log(1-accu)))
}
negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodlin=function(kappaall,kappaloc){
accuvec=exp2acculin(kappaall,kappaloc)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

######################
#linear mutual information no neighborhood

exp2acculinnoneigh=function(kappaall){
accuvec=vector('numeric')
for(i in 1:4){
cost=function(accu){
return((kappaall)*(accu*log(accu)+(1-accu)*log(1-accu)))
}
negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodlinnoneigh=function(kappaall){
accuvec=exp2acculinnoneigh(kappaall)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}


################################
#General Entropy Model

exp2accugen=function(kappaall,kappaloc,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
if(rho==1){
totcost=(kappaall+kappaloc)*(accu*log(accu)+(1-accu)*log(1-accu))
}
if(rho==2){
totcost=(kappaall+kappaloc)*(-1/2*(log(accu)+log(1-accu))-2*log(2))
}
if(rho!=1&rho!=2){
totcost=(kappaall+kappaloc)*(2^(1-rho)/((rho-1)*(rho-2))*(accu*accu^(1-rho)+(1-accu)*(1-accu)^(1-rho))-1/((rho-1)*(rho-2))-log(2))
}
return(totcost)
}


negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodgen=function(kappaall,kappaloc,rho){
accuvec=exp2accugen(kappaall,kappaloc,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

####################################
#Power mutual information

exp2accupow=function(kappaall,kappaloc,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
totcost=(accu*log(accu)+(1-accu)*log(1-accu)-2*0.5*log(0.5))
if(rho>0){
return(totcost^rho)
}else{
return(-1*totcost^rho)
}
}

negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-(kappaall+kappaloc)*cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

exp2likelihoodpow=function(kappaall,kappaloc,rho){
accuvec=exp2accupow(kappaall,kappaloc,rho)
like=exp2multinomconst+sum(exp2totcor*log(accuvec)+exp2totinc*log(1-accuvec))
return(-like)
}

######################################################################################################################################################################################


#import data
symmetrydata=read.dta("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Expt1.dta")
symmetryframe=data.frame(id=symmetrydata$var1,user_id=symmetrydata$var2,qn=symmetrydata$var3,qid=symmetrydata$var4,order=symmetrydata$var5,bid=symmetrydata$var6,chosen_act=symmetrydata$var9, true_state=symmetrydata$var10)

#reward value
r=10 


#data frame for letters
symframeletters=subset(symmetryframe,user_id<1700)
symframeletters$correct=(symframeletters$true_state<=13 & symframeletters$chosen_act==1)+(symframeletters$true_state>13 & symframeletters$chosen_act==2)
unique(symframeletters$qid)

#data frame for balls
symframeballs=subset(symmetryframe,user_id>1700)
symframeballs$true_state=symframeballs$true_state+6*(symframeballs$qid==3)+4*(symframeballs$qid==4)+2*(symframeballs$qid==5)
symframeballs$correct=(symframeballs$true_state<=10 & symframeballs$chosen_act==1)+(symframeballs$true_state>10 & symframeballs$chosen_act==2)
unique(symframeballs$qid)

#Variables to be stored
acculin=matrix(rep(0,80),nrow=4,ncol=20)
accupow=matrix(rep(0,80),nrow=4,ncol=20)
accugen=matrix(rep(0,80),nrow=4,ncol=20)
accunoneigh=matrix(rep(0,80),nrow=4,ncol=20)
dataaccu=matrix(rep(0,80),nrow=4,ncol=20)

#Qids for each treatment type
lettersqid=c(3,5,7,8)
ballsqid=c(3,4,5,6)

#number of states used by different treatments
Nvec=c(8,12,16,20)

#States for different treatments
states=list()
states[[1]]=c(7:14)
states[[2]]=c(5:16)
states[[3]]=c(3:18)
states[[4]]=c(1:20)

correct=list()
incorrect=list()

#{*} Make correct vecs (needs changing between treatment types)
for(i in 1:4){
correct[[i]]=vector('numeric')
incorrect[[i]]=vector('numeric')
for(j in 1:Nvec[i]){
correct[[i]][j]=sum(symframeballs$correct*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
incorrect[[i]][j]=sum((1-symframeballs$correct)*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
}
}



exp4multinomlist=list()
exp4multinomvec=vector("numeric")
for(i in 1:4){
exp4multinomlist[[i]]=vector("numeric")
for(j in 1:Nvec[i]){
exp4multinomlist[[i]][j]=multinom(c(correct[[i]][j],incorrect[[i]][j]))
}
exp4multinomvec[i]=sum((exp4multinomlist[[i]]))
}

exp4multinomconst=sum(exp4multinomvec)

######################
#Linear Mutual Information Costs

exp4accuracylin=function(kappaall,kappaloc,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))
return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodlin=function(kappaall,kappaloc){
accuvec1=exp4accuracylin(kappaall,kappaloc,1)
accuvec2=exp4accuracylin(kappaall,kappaloc,2)
accuvec3=exp4accuracylin(kappaall,kappaloc,3)
accuvec4=exp4accuracylin(kappaall,kappaloc,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
#################################################

#power costs

exp4accuracypow=function(kappaall,kappaloc,rho,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))-N*1/N*log(1/N)
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=1/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))-2*0.5*log(0.5)
}
neighborcost[N/2]=(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))-2*0.5*log(0.5)
return(kappaall*overallcost^rho+kappaloc*sum(neighborcost^rho))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2)-c(1:(N/2))/200,negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodpow=function(kappaall,kappaloc,rho){
accuvec1=exp4accuracypow(kappaall,kappaloc,rho,1)
accuvec2=exp4accuracypow(kappaall,kappaloc,rho,2)
accuvec3=exp4accuracypow(kappaall,kappaloc,rho,3)
accuvec4=exp4accuracypow(kappaall,kappaloc,rho,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

#################################################
#Gen entropy approach
exp4accuracygen=function(kappaall,kappaloc,q,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
neighborcost=vector('numeric')

if(q==1){

overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}else{

if(q==2){
overallcost=kappaall*(-1/N*sum(log(accu/(N/2))+log((1-accu)/(N/2)))-2*log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(-1/2*((accu[j]+accu[j+1])/2*(log(accu[j]/(accu[j]+accu[j+1]))+log(accu[j+1]/(accu[j]+accu[j+1])))+(2-accu[j]-accu[j+1])/2*(log((1-accu[j])/(2-accu[j]-accu[j+1]))+log((1-accu[j+1])/(2-accu[j]-accu[j+1]) ) )) -2*log(2))
}
neighborcost[N/2]=kappaloc*(-1/2*(log(accu[N/2])+log(1-accu[N/2]))-2*log(2))


}else{
overallcost=kappaall*(N^(1-q)/((q-2)*(q-1))*(sum((accu/(N/2))^(2-q)+((1-accu)/(N/2))^(2-q)))-1/((q-2)*(q-1))-log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[j]+accu[j+1])/2*((accu[j]/(accu[j]+accu[j+1]))^(2-q)+(accu[j+1]/(accu[j]+accu[j+1]))^(2-q))+(2-accu[j]-accu[j+1])/2*(((1-accu[j])/(2-accu[j]-accu[j+1]))^(2-q)+((1-accu[j+1])/(2-accu[j]-accu[j+1]))^(2-q)))-1/((q-2)*(q-1))-log(2))

}
neighborcost[N/2]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[N/2])^(2-q)+(1-accu[N/2])^(2-q))-1/((q-2)*(q-1))-log(2))

}

}

return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

exp4likelihoodgen=function(kappaall,kappaloc,q){
accuvec1=exp4accuracygen(kappaall,kappaloc,q,1)
accuvec2=exp4accuracygen(kappaall,kappaloc,q,2)
accuvec3=exp4accuracygen(kappaall,kappaloc,q,3)
accuvec4=exp4accuracygen(kappaall,kappaloc,q,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

######################
#Linear Mutual Information Costs no neighborhoods

exp4accuracynoneigh=function(kappaall,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Information cost function
cost=function(accu){
overallcost=kappaall*(accu*log(accu/(N/2))+(1-accu)*log((1-accu)/(N/2)))
return(overallcost)
}

#Optimization target
negpayoff=function(accuscalar){
return(-1*(r*sum(rep(accuscalar,N/2)/(N/2))-cost(accuscalar)))
}

#Find optimal accuracies
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec=rep(optimout$minimum,N/2)
return(accuvec)
}

exp4likelihoodnoneigh=function(kappaall){
accuvec1=exp4accuracynoneigh(kappaall,1)
accuvec2=exp4accuracynoneigh(kappaall,2)
accuvec3=exp4accuracynoneigh(kappaall,3)
accuvec4=exp4accuracynoneigh(kappaall,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=exp4multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
######################################################################################################################################################################################

likelihoodlin=function(kappaall,kappaloc){
return(exp2likelihoodlin(kappaall,kappaloc)+exp4likelihoodlin(kappaall,kappaloc))
}

modellin=mle(likelihoodlin,start=list(kappaall=3,kappaloc=24),method="L-BFGS-B",lower=c(0.001,0),upper=c(100,100))


likelihoodpow=function(kappaall,kappaloc,rho){
return(exp2likelihoodpow(kappaall,kappaloc,rho)+exp4likelihoodpow(kappaall,kappaloc,rho))
}

modelpow=mle(likelihoodpow,start=list(kappaall=14.1,kappaloc=88.9,rho=2.01),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(1000,1000,1000))

likelihoodgen=function(kappaall,kappaloc,rho){
return(exp2likelihoodgen(kappaall,kappaloc,rho)+exp4likelihoodgen(kappaall,kappaloc,rho))
}

modelgen=mle(likelihoodgen,start=list(kappaall=0.02,kappaloc=1.8,rho=8.4),method="L-BFGS-B",lower=c(0.001,0,0.001),upper=c(100,100,100))

likelihoodnoneigh=function(kappaall){
return(exp2likelihoodlinnoneigh(kappaall)+exp4likelihoodnoneigh(kappaall))
}
modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=4),method="L-BFGS-B",lower=c(0.001),upper=c(100))

#Some LR tests
noneightest=lr.test(-logLik(modelnoneigh),-logLik(modellin),alpha=0.05,df=1)
gentest=lr.test(-logLik(modellin),-logLik(modelgen),alpha=0.05,df=1)
powtest=lr.test(-logLik(modellin),-logLik(modelpow),alpha=0.05,df=1)

#Fits

kappaalllin=coef(modellin)[1]
kappaloclin=coef(modellin)[2]
for(expi in 1:4){
N=Nvec[expi]
accuracy=exp4accuracylin(coef(modellin)[1],coef(modellin)[2],expi)
acculin[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}



kappaallpow=coef(modelpow)[1]
kappalocpow=coef(modelpow)[2]
rhopow=coef(modelpow)[3]
for(expi in 1:4){
N=Nvec[expi]
accuracy=exp4accuracypow(coef(modelpow)[1],coef(modelpow)[2],coef(modelpow)[3],expi)
accupow[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}

kappaallgen=coef(modelgen)[1]
kappalocgen=coef(modelgen)[2]
rhogen=coef(modelgen)[3]

for(expi in 1:4){
N=Nvec[expi]
accuracy=exp4accuracygen(coef(modelgen)[1],coef(modelgen)[2],coef(modelgen)[3],expi)
accugen[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}

kappaallnoneigh=coef(modelnoneigh)[1]

for(expi in 1:4){
N=Nvec[expi]
accuracy=exp4accuracynoneigh(coef(modelnoneigh)[1],expi)
accunoneigh[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}


###################################################################
#Create the actual data with bootstrap standard errors
#Store upper and lower bounds of 95% CI
dataaccuupper=matrix(rep(0,80),nrow=4,ncol=20)
dataacculower=matrix(rep(0,80),nrow=4,ncol=20)

#number of bootstrap iterations
rep=1000

for(expi in 1:4){
#Setup for speicifc treatment
playframe=subset(symframeballs,qid=ballsqid[expi])
N=Nvec[expi]

#Define Players
players=unique(playframe$user_id)
nplay=length(players)

#Set up value storage
acculist=list()
for(i in 1:N){
acculist[[i]]=vector("numeric")
}

#Draws
drawvec=sample(c(1:nplay),nplay*rep,replace=TRUE)

for(i in 1:rep){
#define bootframe
bootframe=data.frame(id=vector("numeric"), user_id=vector("numeric"), qn=vector("numeric"), order=vector("numeric"), bid=vector("numeric"),chosen_act=vector("numeric"), true_state=vector("numeric"), correct=vector("numeric"))

#fill bootframe
for(j in 1:nplay){
bootframe=rbind(bootframe, subset(playframe, user_id==players[drawvec[(i-1)*nplay+j]] ) )
}

#caluculate accuracies
for(j in 1:N){
bootsubframe=subset(bootframe,true_state==states[[expi]][j])
acculist[[j]][i]=sum(bootsubframe$correct)/length(bootsubframe$correct)
}

}

accuracyupper=vector("numeric")
accuracylower=vector("numeric")

#sort vectors
for(i in 1:N){
acculist[[i]]=sort(acculist[[i]])
accuracyupper[i]=acculist[[i]][975]
accuracylower[i]=acculist[[i]][25]
}

#calculate
accuracy=correct[[expi]]/(correct[[expi]]+incorrect[[expi]])

#fill in accuvec
dataaccu[expi,]=c(rep(0,(20-N)/2),accuracy,rep(0,(20-N)/2))
dataaccuupper[expi,]=c(rep(0,(20-N)/2),accuracyupper,rep(0,(20-N)/2))
dataacculower[expi,]=c(rep(0,(20-N)/2),accuracylower,rep(0,(20-N)/2))
}

#Table A4.3************************************************************************************************
#some output tables
dp11frame=data.frame(state=states[[1]],data=dataaccu[1,states[[1]]],datahigh=dataaccuupper[1,states[[1]]],datalow=dataacculower[1,states[[1]]],mutual_info=acculin[1,states[[1]]],power_info=accupow[1,states[[1]]],gen_ent=accugen[1,states[[1]]],no_neigh=accunoneigh[1,states[[1]]])
stargazer(dp11frame,summary=FALSE)
write.table(dp11frame, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/ballsfits_combo.csv", sep = ",", col.names = T, append = F)

dp12frame=data.frame(state=states[[2]],data=dataaccu[2,states[[2]]],datahigh=dataaccuupper[2,states[[2]]],datalow=dataacculower[2,states[[2]]],mutual_info=acculin[2,states[[2]]],power_info=accupow[2,states[[2]]],gen_ent=accugen[2,states[[2]]],no_neigh=accunoneigh[2,states[[2]]])
stargazer(dp12frame,summary=FALSE)
write.table(dp12frame, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/ballsfits_combo.csv", sep = ",", col.names = T, append = T)

dp13frame=data.frame(state=states[[3]],data=dataaccu[3,states[[3]]],datahigh=dataaccuupper[3,states[[3]]],datalow=dataacculower[3,states[[3]]],mutual_info=acculin[3,states[[3]]],power_info=accupow[3,states[[3]]],gen_ent=accugen[3,states[[3]]],no_neigh=accunoneigh[3,states[[3]]])
stargazer(dp13frame,summary=FALSE)
write.table(dp13frame, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/ballsfits_combo.csv", sep = ",", col.names = T, append = T)

dp14frame=data.frame(state=states[[4]],data=dataaccu[4,states[[4]]],datahigh=dataaccuupper[4,states[[4]]],datalow=dataacculower[4,states[[4]]],mutual_info=acculin[4,states[[4]]],power_info=accupow[4,states[[4]]],gen_ent=accugen[4,states[[4]]],no_neigh=accunoneigh[4,states[[4]]])
stargazer(dp14frame,summary=FALSE)
write.table(dp14frame, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/ballsfits_combo.csv", sep = ",", col.names = T, append = T)
#*******************************************************************************************************


#Table 4.2****************************************************************************
fitnoneigh=exp2acculinnoneigh(kappaallnoneigh)
fitnoneigh
fitlin=exp2acculin(kappaalllin,kappaloclin)
fitlin
fitgen=exp2accugen(kappaallgen,kappalocgen,rhogen)
fitgen
fitpow=exp2accupow(kappaallpow,kappalocpow,rhopow)
fitpow
#*************************************************************************************
